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Abstract  _  fe  l\  . 


Suppose  that  two  stationary  time  series  satisfy  the  linear 
relationship  y(m  )  ■  J  m«gh(mr  )x((n-m)r  ). ^Estimating  the  distributed  lag 


coefficients  {h(mr)}  from  a  sample  of  the  two  processes  when  (x(nr))  and 

V  rTTWO  ) 

{ y(nr )}  are  measured  with  error  is  a  statistical  problem  that  is 
frequently  encountered  in  physical  science,  engineering,  and  social 
science  applications.  In  the  engineering  and  science  literature  the 
distributed  lags  are  called  the  impulse  response  weights  of  a  causal 
linear  filter.  A  least  squares  fit  of  the  model  gives  biased  estimates 
of  the  coefficients  for  this  time  series  version  of  the  errors-in- 
variables  problem.  This  paper  presents  approximately  unbiased  estimators 
of  a  scalar  multiple  of  the  coefficients.^The  large  sample  variances  of 
these  estimators  are  identical,  and  are  of  the  order  OCN”1)  where  N  is 
the  sample  size.  The  accuracy  and  practicality  of  the  estimating 
procedure  is  illustrated  by  some  simulation  results. 
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Introduction 

Suppose  that  two  discrete-time  stationary  random  processes  satisfy 

00 

the  linear  relationship  y(tn)  »  £  h(tm)x(tn  -  tm)  where  tn  -  nr . 

m*-*> 

Estimating  the  sequence  {h(tm)}  from  a  sample  of  the  two  processes  is  a 
statistical  problem  that  is  frequently  encountered  in  physical  science, 
engineering,  and  social  science  applications.  Often  (x(tn)}  and  {y(tn)} 
are  measured  with  error,  and  are  of  the  form  x(tn)  -  x(tn)  +  u(tn)  and 

•m 

y(tn)  »  y(tn)  +e(tn)  where  (e(tn)}  and  {u(tn)>  are  noise  processes. 

Then  a  least  squares  fit  of  the  model  yields  biased  estimates  of  the 
h(tn).  These  biases  are  examples  of  the  error s-in-varlables  biasing 
problem  for  a  linear  model  (Kmenta,  Chapter  9,  1971).  The  bias  of  the 

A 

least  squares  estimator  h(tn)  varies  with  n  unless  the  covariance  of 

A 

(u(tn)}  is  identical  to  that  of  (x(tn)}  .  The  pattern  of  {h(tn»  will  be 
completely  different  from  the  pattern  of  the  true  h(tn)  if  the  errors  are 
sizable  as  compared  with  the  true  series  (Clay  and  Hinich,  1981). 

Assuming  some  mild  restrictions  on  (h(tn)}  and  the  additive  errors. 
Clay  qnd  Hinich  and  Hinich  (1982)  develop  consistent  estimators  of  a 
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This  paper  considers  the  discrete-time  linear  relationship  rather 
than  the  continuous-time  model  used  in  the  cited  previous  work.  A 
Hilbert  transform  relationship  for  sequences  is  used  to  derive  consistent 
estimators  of  the  h(tn)  from  a  consistent  estimator  of  the  cross  spectrum 

«*  m 

between  (x(tn)}  and  {y(tn)} .  The  approach  uses  the  frequency  domain 
Interpretation  of  a  linear  filter  (Brilllnger,  Section  2.7,  1975).  The 
basic  terminology  necessary  to  understand  the  method  is  reviewed  in  the 
next  section. 

1.  The  Discrete  Hilbert  Transform 

In  filter  theory  terminology »  (x(t)}  is  the  input  to  a  linear  filter 
whose  Impulse  response  is  {h(tn»  .  The  filter  is  called  causal  if  h(tn) 

■  0  for  tn  <  0,  and  thus  y(tn)  depends  only  on  x(sn)  for  sn  <  tn.  In 
economics ,  a  causal  linear  model  is  called  a  distributed  lag  with  lag 
parameters  {h(tn)J  . 

00 

A  filter  is  stable  if  £  n-ol^^n)  |<  00  •  The  frequency  response  of  a 
stable  filter  is  characterized  by  its  transfer  function 
00 

H(f )  -  l  h(tn)exp(-12iT  ftn).  (1.1) 

n-0 

The  inverse  relationship  is 
1/t 

h(tn)-TJ  H(f)exp(i2irftn)df.  (1.2) 

0 

Since  ■  nx  ,  H(f  +  1/t  )  ■  H(f)  for  all  f.  For  real  {h(tn)>  ,  H(-f)  - 
H*(f)  where  the  star  denotes  complex  conjugate,  and  thus  the  phase 
response  ♦(f)  ■  tan“l{ ImH(f )/ReH(f ) ]  is  an  odd  function  with  <fr(0)  ■  0  for 


the  range  -it  <  $  <  x  .  The  filter's  gain  | H(f )  |  is  an  even  function. 


The  Hilbert  transform  relates  ln|H(f)|  to  $(f)  if  the  causal  filter 

00 

is  minimum  phase-lag,  i.e.  if  H(z)  h(tn)zn  ^as  no  zeros  on  or 

inside  the  unit  circle  |z|  »  1  on  the  complex  plane.  The  phase-lag  -^(f) 
of  a  finite  minimum  phase-lag  filter  is  less  than  any  other  filter  with 
the  same  gain  function  (Zadeh  and  Desoer,  Section  9.7,  1963).  Normalize 
the  time  unit  so  that  t  ■  1.  The  log  gain  is  related  to  its  phase  by  the 
Hilbert  transform 

log|H(f)|  -  J1  <J>(g)cot[ir  (f-g)ldg  +  c,  (1.3) 

0 

where  c  is  a  scaling  constant  (Gold  and  Rader,  p.  248,  1969).  Note  that 
the  integrand  has  a  singularity  only  at  g  ■  f  where  the  principal  value 
of  the  integral  exists. 

Suppose  that  the  gain  of  a  minimum  phase-lag  filter  is  normalized  by 
setting  c  -  0.  Then  it  is  proven  in  the  Appendix  that  h(0)  ■  ±1. 

It  is  clear  from  (1.3)  that  accurate  estimates  of  the  phases  at  a 
dense  frequency  grid  on  the  unit  interval  yield  accurate  estimates  of  a 
scalar  multiple  of  the  gain  at  these  frequencies  using  a  finite  sum 
approximation  to  the  integral.  The  estimated  gain  combined  with  the 
estimated  phases  gives  an  estimate  of  the  transfer  function  which  can  be 
inverted  to  yield  an  estimate  of  the  impulse  response  (distributed  lag) 
coefficients.  If  the  estimated  impulse  response  is  minimum  phase,  then 
it  follows  from  the  theorem  in  the  Appendix  that  the  estimate  of  h(0)  is 
nearly  one  when  h(0)  ■  1  for  the  normalization  c  «  0,  provided  the  errors 
are  sufficiently  small  due  to  the  use  of  a  large  sample  or  because  the 
signal-to-noise  ratio  in  the  data  signals  are  large. 


Accurate  estimates  of  the  phases  can  easily  be  derived  from  a 
consistent  estimator  of  the  cross  spectrum.  Assume  that  the 
autocovariance  of  {x(tQ)}  is  absolutely  summable.  Then  the  spectrum 
S^f)  of  {x(tn}}  exists,  and  more  importantly,  the  cross  spectrum  is  of 
the  form 

Sxy(f)  -  HCfJSxCf)  (2.1) 

(Jenkins  and  Watts,  Sect.  8.4.2,  1968).  If  S^f)  >  0,  it  follows  from 
(2.1)  that 

4(f)  -  arctan[ImSxy(f)/ReSxy(f)].  (2.2) 

Thus  a  consistent  estimator  of  Sxy(f)  yields  a  consistent  estimator  of 
4(f)  using  (2.2). 

Recall  that  the  sample  consists  of  simultaneous  observations  of  the 
•  00 

noisy  processes  {x(tn)  -  x(tn)  +  u(tn)>  and  (y(tn)  -  y(tn)  +  e(tn)> .  The 

following  assumptions  are  crucial  for  the  method:  1)  {u(tn)>  and  {e(tn)> 

are  uncorrelated,  and  2)  they  are  uncorrelated  with  the  true  processes 

{*(tn)}  and  {y(tn)J  .  Given  these  assumptions,  the  cross  covariance 

functions  of  the  observed  and  true  processes  are  Identical,  l.e.  c- - ( t )  - 

xy 

cxy(t)  for  all  t.  Thus  the  cross  spectrum  of  the  observed  processes  is 

SxyU). 

There  are  several  related  approaches  to  estimating  the  cross 
spectrum  of  the  observed  processes.  Consistent  estimators  and  their 
asymptotic  properties  are  presented  by  Anderson  (Chapter  9,  1971), 
Brilllnger  (Chapter  5,  1975),  and  Fuller  (Chapter  7,  1976). 
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One  simple  method  uses  the  discrete  Fourier  transforms  of  a  sample 
«* 

(x(tn),  y(tn)|n  -  0, 1, . . . ,N-1}  .  For  frequencies  fk  -  k/N,  these 
transforms  are: 


*  N-l  ~ 

X(fk)  "  I  x(tn)exp(-12rnk/N) 
n-0 


and 


(2.3) 


~  N-l  - 

Y<fk)  ■  I  y(tn)exp(-i2rnk/N). 
n*0 


The  cross  spectrum  estimator  at  frequency  0  <  fQ  <  1/2  is  the  average^ 


-  I  (M-l)/2 

S~~(f0)  “  —  l  X(fQ  +  fk)Y*(fc  +  fk),  (2.4) 

xy  MN  k— (M-D/2 


where  M«N.  From  Theorems  7.3.1  and  7.3.2  in  Brillinger  (1975),  the 

A 

expected  value  of  S-~(f0)  is  S»-(f0)  +  0(M/N)  and  Its  variance  is  0(M-1). 

xy  xy 

Thus,  the  mean  square  error  is  small  if  and  only  if  M  is  large  and  M/N 
is  small.  To  parameterize  the  tradeoff  between  bias  and  variance,  let 
M  *  N*  where  0  <  a  <  1.  Good  results  results  for  the  mean  square  error  of 

A 

the  phase  of  S-~(f0)  have  been  obtained  in  a  variety  of  simulated 
Xy 

distributed  lag  models  that  have  tried  using  ft<Nl/2  and  N*500. 

The  estimate  of  phase  is  given  by 


$(f0)  -  arctanf  Im  S--(fQ)/Re  S--  (f 0)  ] . 

xy  xy 


(2.5) 
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Expanding  (2.5)  in  a  Taylor  series  about  S--(f0),  it  follows  that 
*  xy 

4>(f0)  +  e  where  the  expected  value  of  the  error  e  is  of  order  M/N. 

A  A 

For  1/2  <  fQ  <  1,  <f>(£0)  m  -$(l-fQ).  The  variance  of  the  phase  error  is 

a^(fQ)  -  (2M)-1  W ”2Cf0)  -  1]  +  0(M-2),  (2.6) 

where 

Y(f0>  “  |S--(f0)|/[S-(f0)S-(f0)]l/2  (2.7) 

xy  x  y 

denotes  the  coherence  between  the  observed  signals. 

A 

Since  M  »  IP  ,  it  follows  from  (2.6)  and  the  above  that  $(f0)  a 

A 

consistent  estimator  of  <p  (f0)  as  if  y(fQ)  >  O.2  Also,  <|>(fj)  and 

A 

<fr(fk)  are  approximately  uncorrelated  if  | j  —k. |  >  M  because  the  discrete 
Fourier  coefficients  have  cross  correlations  of  order  N“1  if  the  two  time 
series  have  well  behaved  cumulants  (Brillinger,  Theorem  4.4.2,  1975). 

A 

If  so,  <fr(^k(N))  is  asymptotically  Gaussian  if  k(N)/N*-f0  as  N*-«,  eg.  k(N)  - 
[fQN]  where  the  brackets  denote  the  integer  part  of  the  number. 

Expression  (2.6)  can  be  rewritten  in  terms  of  the  noise-to-signal 
ratios  for  the  observed  processes.  These  ratios  are  denoted  rx(f)  - 
Su(f)/Sx(f)  and  ry(f)  -  S^(f)/Sy(f),  where  Su(f)  and  Sfe(f)  are  the 
spectra  of  the  respective  additive  noises.  Then  from  (2.7),  y2(f)  » 
[(l+rx(f))(l+ry(f))l-l.  Thus 

2Mo*(f)  -  rx(f)  +  ry(f)  +  rx(f)ry(f)  +  0(M-1), 


(2.8) 


which  implies  that  is  approximately  linear  in  the  noise-to-signal 

2 

ratios  when  they  are  small.  The  coherence  becomes  small,  and  o ^  large, 
when  rx  and  ry  are  big  (low  signal-to-noise  ratios). 

3.  Estimating  the  Gain 

A  finite  sum  approximation  to  the  integral  in  (1.3)  for  f  -  k/N  yields 
N-l 

log|H(f)|  -  N-l  l  ♦(gj)cot[ir(fk  -  gj)]  +  0(N-1)  +  c,  (3.1) 
j-0 
j*k 

where  the  sum  omits  j*k  to  avoid  the  singularity  at  j«k.  Substituting  the 

A 

estimates  <t>(gj)  for  phases  in  (3.1)  yields  a  simple  estimate  of 

A 

log|H(fjc)|  +  c  which  will  be  denoted  esttlogl^f^)  |  ] .  Since  the  $(gj)  are 
Gaussian  and  consistent,  it  follows  from  (3.1)  that  esttlogl^f^^jjj)  |  ]  is 
a  consistent  Gaussian  estimator  of  log|H(f0)|  as  N»-<».  The  large  sample 
variance  of  this  estimator  will  now  be  derived. 

Theorem  1.  The  large  sample  variance  of  esttlogl^f^jj))  |  ]  is 

b  2 

o2(*k(N))  -  *”2  I  ®"^(J»(fk(N)  +  SmM>»  (3.2) 

m—a+1 
n*  0 


where  a  -  [f0N/M]  and  b  -  [(l-f0)N/M],  and  g,^  -  mM/N.  Thus  the 
large  sample  standard  deviation  is  of  order  if  ct-l/2  (a  resolution 

bandwidth  of  N"l/2). 

i-f 

Proof:  The  Integral  in  (1.3)  can  be  rewritten  as  -  J  $  (f+y)cot(iry)dy. 
A  finite  sum  approximation  of  this  Integral  that  converges  as 


when  is 
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b 

-(M/N)  l  <|>  (ffc  +  ggiM)cot(ir  gmM)  •  (3.3) 

m*-a+l 
m*  0 


Thus  the  large  sample  variance  of  (3.3)  with  $  replacing  $  Is 
asymptotically  equal  to  the  large  sample  variance  of  est[log|H(fk)  j  ] .  The 

A  A 

correlation  between  $ (mM/N)  and  $ (nM/N)  is  0(N“1)  for  n*n,  and  thus  It 
follows  from  (2.6)  that  the  variance  of  (3.3)  for  large  N  is 
b  2 

(M/N)2  £  <>$(flc  +  gmM^o^^SmM)  •  (3.4) 

m— a+1 
m*  0 

Since  cot^  ■  sin“2x  -  1  and  N-2  sin“2(ir  mM/N)  ■  (ttoM)-2  +  0(N“2),  for 

b  2 

large  N  expression  (3.4)  becomes  ir“2  £  m**fy  (f k  +  g,^) ,  the  desired 

m— a+1 
n*0 

large  sample  variance  o2(fic). 

For  example,  suppose  that  the  coherence  is  constant  over  frequencies. 
00 

Since  £  m"2  -ir2/6,  then  the  large  sample  variance  is  (y“2  -  1)/6M. 

m-1 

When  y(f)  is  slowly  varying,  a  useful  approximation  of  <j2(fk)  is 

h"2(fk)  -  1]/6M. 

It  will  now  be  shown  that  the  estimates  are  asymptotically 
uncorrelated. 


Theorem  2.  The  large  sample  covariance  between  est[log|H(f j (Nj) | ]  and 
est[log|H(fk(N))|  J  for  j(N)  -  Cf  J.N]  and  k(N)  -  [f2N]  is  0(1/N|f  2-f  i| ). 


Proof:  Once  again  the  asymptotic  Independence  of  the  phase  estimators 
will  be  used.  From  (3.3)  and  (3.4),  the  large  sample  covariance  is 


approximated  by 
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2 

(M/N)2  £  +  gnMJcot^gmM)  cotMgnji  -  t<i  +  f  i>  J  (3.5) 

m 

if  fj  <  f 2»  The  sum  excludes  Che  singularities  at  m  -  0,  and  m  - 
(f2“fi)N/M  if  it  is  an  integer  in  the  set  {-  [f^N/M]  <  m  <  [ ( 1— f 2) N/M]} . 
The  largest  values  of  the  summand  are  at  m  -  ±  1  and  m  ■  [(f 2~f i)N/M]±  1, 
where  they  are  approximately  (H/N)2m“1(N/M)/it  2(f 2-f  j)  ■  [u  2N(f2*-f  i)]”1. 

4.  Comparison  with  Other  Methods 

The  numerical  approximations  of  the  Hilbert  transform  described  in 
Clay  and  Hinich  and  Hinich  can  yield  very  biased  estimates  in  some  cases 
when  the  sample  size  is  moderate  and  the  bandwidth  is  large.  The  method 
presented  here  usually  gives  less  finite  sample  bias  than  these  other 
methods,  especially  for  discrete-time  filters. 

As  is  to  be  expected  however,  the  large  sample  variance  of 
estlloglHCffc) | ]  is  somewhat  larger  than  that  of  the  other  estimators.  To 
compare  with  the  older  methods ,  suppose  that  the  cross  spectrum  is 
estimated  from  P  successive  samples  of  size  N  (observation  time  Nr )  of 

<*»  m 

(x(tn),  y(tQ)} .  The  sample  cross  spectrum  at  f  is  the  average 

-  "  * 

(p-l,...,P)  of  Che  unsmoothed  terms  Xp(f)Yp(f).  If  the  samples  are 

approximately  uncorrelated,  the  large  sample  variance  of  the  phase 

.  2 

computed  from  this  block  smoothed  cross  spectrum  is  P"1o(^(f)  (Hinich  and 
Clay,  1968).  Reworking  the  proof  of  the  theorem  using  the  estimator 
l(l-f0)Nl  - 

-N-*1  \  ♦  ((k+j  )/N)cot(ir  j/N) ,  its  large  sample  variance  is  0(P“1). 

j— [f0N]+l 

The  large  sample  variance  of  the  other  methods  is  of  order  o(P~l)  as 


5.  Estimating  the  Transfer  Function  and  Impulse  Response 

Recall  that  the  log  gain  is  estimated  up  to  an  additive  scale 

constant.  Thus  exp{ est[log|H(fk) | ]}  is  a  consistent  estimator  of  c|H(fk)| 

Using  the  linear  approximation  of  the  exponential.,  the  variance  of  this 

estimator  is  approximately  |H(fk) 2(fk)  for  large  M.  The  covariance 

*  ,  2 

between  this  estimator  and  <fr(fk)  is  approximately  ir  (fk)  from  (3.3). 
The  transfer  function  is  estimated  by 


H(fk)  -  exp{  est[log| H(fk)  |  ]  +  i  (flc)  h  •  (5.1) 

Since  the  error  in  the  phase  estimate  is  additive,  the  error  in 

A 

est  log|H(fk)|  is  additive  and  thus  the  error  in  H(fk)  is  multiplicative. 

The  complex  variance  of  the  exponent  in  (5.1)  is  defined  to  be  the 

sum  of  the  variances  of  the  real  and  imaginary  parts,  which  is 
2 

<j2(fk)  +  »  (2/3M)  [y  "2(fk)-i J  using  the  simplifying  approximation  for 

<j2(fk)  previously  presented.  Using  the  linear  approximation  of  the 
exponential  function,  for  large  H  the  complex  variance  of  the  estimator  is 

A 

E|H(fk)  -  H(fk)|2  «  (2/3M)|H(fk)|2[Y-2(fk)_i].  (5.2) 

A 

From  Theorem  2  and  the  properties  of  the  $ (fk) ,  the  complex  covariance 

A  A 

between  H(fj)  and  H(fk)  is  0(l/N|j-k|). 

The  estimator  of  the  impulse  response  coefficient  h(tn)  is 


h(tn)  -  N-l  l  H(fk)exp(i2rkn/N).  (5.3) 

k-0 


If  the  estimated  impulse  response  is  minimum  phase,  then  h(0)*± 1+0(N“1) . 
The  error  is  only  due  to  the  errors  from  approximating  the  integrals. 

A 

For  large  M,  the  variance  of  h(tn)  is  approximated  by  the  variance  of 
[N/M]  - 

(M/N)  £  HCfnj^expCiarmM/N) ,  with  an  error  of  0(M/N).  This  variance, 

m*0  a 

with  an  error  of  0(1/NM)  due  to  the  covariances  of  the  HCf,^),  is 

[N/M]  -  [N/M] 

(M/N)2  l  ElHCf^)  -  HCWp  »  (2M/3N2)  £  (HCf^) 

m*0  m*0 


(5.4) 

2  1 

-  ~  /  |H(f)|2[Y-2(f)_ijdf. 

3N  o 


The  largest  errors  in  these  approximations  are  0(M“2)  or  0((M/N)~),  and 
thus  are  0(N~6)  where  0  »  2min(a,l-a).  Note  that  this  variance  is 
Independent  of  n,  so  all  parameter  estimates  have  equal  variance  for  large 
N. 

As  an  example,  suppose  that  coherence  is  constant  across  frequencies. 
1 

Since  /  | H( f ) | 2  df  ■  £ q-q  h2(tn)  by  Parceval's  formula,  assuming 
0 

existence,  then  from  (2.8)  and  (5.4)  the  large  sample  variance  of  h(tn) 
is 


Var  h(tn) 


2 

—  (rx  +  ry  +  rxry) 


00 

l  h2(tn). 
n-0 


(5.5) 


Thus  h(tn)  is  a  precise  estimator  when  N  is  large  and  the  phase  error 
is  small.  Expression  (5.5)  overstates  the  error  in  the  estimates  when  the 


right  hand  side  is  small  since  it  is  derived  without  incorporating  the 

*  N-l  A 

constraint  given  by  h(0)  -  N“l£  ^-0  H(fk)  ■  ±  1  +  0(N-1). 

6.  Data  Analysis 

The  method  was  tested  using  artificially  created  data  for  the  causal 

minimum  phase-lag  filter  {h(0)  -  1,  h(l)  »  3/2,  h(2)  -  1,  h(3)  -  1/2, 

h(4)  ■  1/4,  h(m)  ■  0  for  m  >  5} .  For  each  trial  a  sequence  of  N  -  500 
4 

values  of  y(n)  »  £m„gh(m)x(n-m)  were  computed  using  these  five  lag 
parameters  and  an  input  sequence  generated  by  the  recursion 

x(n)  -  0.5  x(n-l)  +  z(n),  (6.1) 

where  the  z(n)  values  were  generated  by  the  IMSL  subroutine  GGNML  to 
simulate  a  sequence  of  independent  normal  N(0,1)  variates.  In  other 
words,  the  input  to  the  filter  was  a  pseudo-  AR(1)  process  with  parameter 
1/2  and  variance  4/3.  Two  independent  white  noise  sequences  were 
generated  to  produce  noisy  "observations"  x(n)  *  x(n)  +  ou(n)  and 
y(n)  -  y(n)  +oe(n)  where  the  positive  parameter  a  controls  the  standard 
deviation  of  the  errors.  The  pseudo-random  variates  u(n)  and  e(n)  were 
also  N(0, 1)  as  generated  by  the  GGNML  subroutine  using  different  seed 
values.  A  variety  of  seed  values  were  used  to  test  the  method  for  each 
set  of  parameter  values  that  were  tried. 

Table  1  presents  typical  results  for  o  -0.71,  1.00,  1.22,  1.41,  and 
1.58  (o^  ■  0.5,  1,  1.5,  2,  2.5)  using  M  -  39  as  the  smoothing  number.  Th 
resolution  bandwidth  is  then  M/N  -  0.078.  The  same  artificial  data  was 

A 

used  in  each  run.  Table  1  presents  h(t)  for  t  -  0,...,6;  the  root  mean 

A 

square  error  (RMSE)  of  h(t)  for  t  »  0,...,499;  the  square  root  of  the 
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large  sample  variance  (5.4)  with  Che  estimaCed  gain  and  squared  coherence 

A 

in  place  of  Cheir  true  values  (Gst  SD  h(t));  and  Che  average  squared 
coherence  for  0  <  f  <  1/2.  The  coherence  is  large~»-  for  low  frequencies 
since  (x(t)}  is  red  noise  whereas  Che  additive  errors  are  white. 

The  method  produced  estimates  of  the  distributed  lags  that  track  the 
impulse  response  even  for  the  largest  value  of  a  that  was  used.  The 
estimates  of  the  zero  tail  of  the  impulse  response  are  very  accurate, 
especially  for  the  far  right  tail.  For  example,  the  maximum  error  for 
t  >  450  is  0.0004  for  a  2  <  1.5,  0.003  for  a  2  ■  2,  and  0.01  for  o2  ■  2.5. 
This  accuracy  is  greater  than  that  predicted  by  the  asymptotic 
calculations.  When  the  phase  error  is  small,  the  length  of  the  estimated 
impulse  response  is  slightly  larger  than  the  true  length. 

The  estimated  standard  deviation  has  an  upward  bias  due  to  biases  in 
the  estimated  phase  and  gain.  The  large  sample  variance  appears  to  only 
serve  as  an  upper  bound  for  the  sample  size  N  -  500.  The  largest  errors 
are  in  the  estimates  of  h(l)  -  1.5.  These  estimates  are  downward  biased 
by  the  convolution  of  the  true  lags  with  the  transform  of  the  errors  in 
the  estimated  gain.  The  bias  is  a  function  of  M. 

Table  2  presents  results  for  o  -  1.22  using  M  -  9,  19,  29,  39,  49, 
and  59.  The  results  show  how  important  it  is  to  use  different  values  of  M 
when  using  this  method  to  estimate  the  lags.  More  research  is  needed  to 
determine  practical  ways  for  setting  M  as  a  function  of  prior  beliefs 
about  the  lag  pattern,  the  signal-to-noise  ratios,  and  the  shape  of  the 
signal  spectra. 
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Table  3  presents  ordinary  least  squares  estimates  of  the  h(t)  for 
o^  ■  l.S  using  fifteen  lags  and  an  intercept  constrained  to  be  zero.  The 
estimates  are  downward  biased  as  expected.  Although  the  number  of 
non-zero  lags  is  correctly  determined  by  inspection  of  the  t  ratios,  the 
lag  pattern  is  distorted  in  the  OLS  estimates. 

7.  Conclusion 

The  artificial  data  results  indicate  that  the  Hilbert  transform 
method  can  become  a  powerful  tool  for  estimating  distributed  lags  when  the 
signals  are  observed  with  additive  noise.  This  estimation  procedure  can 
give  accurate  estimates  of  the  distributed  lag  coefficients  when  there  is 
considerable  additive  noise  in  the  observed  time  series,  provided  that  the 
sample  size  is  sufficiently  large  to  yield  accurate  phase  estimates.  If 
the  phase  of  the  transfer  function  can  be  estimated  sufficiently  closely, 
then  the  Impulse  response  function  is  determined  by  the  Hilbert  transform 
up  to  a  scalar  multiple.  The  statistical  computational  aspect  of  the 
procedure  is  simple. 
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Footnotes 


If  0  <  f„  <  M/2N,  then  average  the  X*(f0  +  f^JYCfg  +  f^)  for 
0  <  fjt  <  f0  +  M/2N.  A  similar  constriction  of  bandwidth  holds  when 
1/2-M/2N  <  fQ  <  1/2. 


A  consistent  estimator  ofy(fQ)  is  |S-»  (f0)  |/ [S-  (fQ)S-  (fQ)  ] 

.  xy  x  y 

where  S-(fQ)  and  S>(fQ)  are  consistent  estimators  of  the  two  spectra 
x  y 

A  prudent  policy  is  to  not  estimate  the  phase  at  fr  a  r/N  if  its 
associated  estimated  coherence  is  below  a  preset  threshold.  The 
frequencies  where  the  coherence  is  too  low  are  skipped  in  the  finite 
sum  approximation  of  the  integral  in  (1.2). 
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Table  3 

OLS  Estimates  for  Fifteen  Lags  and  a 
Zero  Intercept 


qi  m  1<5 

Variable  Estimate  Standard  Error  t  Ratio 


h(0) 

0.53 

0.07 

7.90 

h(l) 

0.69 

M 

9.95 

h(2) 

0.6A 

9.25 

h(3) 

0.A6 

6.65 

h(A) 

0.30 

A.2A 

h(5) 

0.12 

1.78 

h(6) 

0.0A 

0.62 

h(7) 

0.0A 

0.5A 

h(8) 

0.07 

0.9A 

h(9) 

0.03 

0.A0 

h(  10) 

-0.05 

-0.68 

h(  11) 

0.01 

0.15 

h<  12) 

-0.05 

-0.67 

h(13) 

-0.06 

-0.92 

h(lA) 

-0.03 

-O.AA 

h(15) 

0.00 

0.0A 

R2  -  0.5A  with  and  without  Intercept 
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APPENDIX 

m 

If  It-0lh<O|  <  • ,  the  analytic  continuation  of  the  transfer 
function  H(f)  inside  the  unit  circle.|z|  -  1  in  the  complex  z-plane  is 

m 

H<.(z)  “  I 0°  unit  circle  z  ■  exp(-12irf),  H^z)  ■  H(f). 

For  z  -  0,  1^(0)  -  h(0). 

If  H(f)  is  minimum  phase,  H^z)  has  no  zeros  in  |z|  <  1  and  the 
Hilbert  transforms  of  phase  and  gain  are: 

♦(f)  -  /1log|H(g)|cot  x (g-f)dg  (Al) 

0 

and 

log|H(f)|  -  c  +  /  Vg)cot  x  (f-g)dg,  (A2) 

0 

where  c  is  an  arbitrary  constant.  Set  c  -  0. 

The  following  lemmas  are  used  In  the  proof  of  the  main  result. 

,1 

Lemma  _1.  J  cot  x  (y-*)dy  ■  0  for  every  x. 

~  0 

Proof:  Consider  the  minimum  phase  filter  defined  by  {h(0)  ■  2,  h(t)  -  0 
for  t  *  0} .  This  filter’s  transfer  function  is  the  constant  2,  and  thus 
its  phase  is  zero  for  all  f.  From  (Al), 

0  ■/  log|  2|cot  x  (y-x)dy, 

0 

and  the  result  then  follows. 
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Lemma  2.  /  log|H(f)|df  ■  0  for  a  minimum  phase  filter. 

0 

Proof:  Reversing  the  order  of  integration  when  Integrating  both  sides  of 
(A2) , 

/  log|H(f)|df  -/  [/  *cot  ir  (f-g)df]<fr  (g)dg  -  0 
0  0  0 

from  Lemma  1. 

Theorem.  h(0)  -  ±  1. 

Proof:  Since  H(z)  has  no  zeros  or  poles  in  |z|  <  1,  by  Jensen's  theorem 
(p.  125,  Titchmarsh,  1952) 

log |H(0)  |  - /^oglHCfJldf. 

0 


From  Lemma  2,  log[H(0)(  -  0  and  thus  log(h(0)l  »  0. 
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